Insurance
mydata <- read.table("Insurance.csv", header=TRUE, sep=",", dec=".")
mydata <- mydata[-6]
mydata$ID <- seq(1, nrow(mydata))
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmydata <- mydata %>%
mutate(Gender = replace(Gender, Gender == "female", "F")) %>%
mutate(Gender = replace(Gender, Gender == "male", "M")) %>%
mutate(Smoker = replace(Smoker, Smoker == "yes", "Yes")) %>%
mutate(Smoker = replace(Smoker, Smoker == "no", "No"))
set.seed(10)
mydata <- mydata[c(7, 6, 1, 2, 3, 4, 5)] [sample(1:1338, size = 140), ]
mydata$Cost <- round(mydata$Cost, 0)
mydata$BMI <- round(mydata$BMI, 1)
rownames(mydata) <- seq(1, nrow(mydata))
head(mydata) ## ID Cost Age Gender BMI Kids Smoker
## 1 491 1749 19 F 32.9 0 No
## 2 368 8017 42 F 25.0 2 No
## 3 439 12593 52 F 46.8 5 No
## 4 344 13982 63 M 36.8 0 No
## 5 1295 11931 58 M 25.2 0 No
## 6 143 18972 34 M 25.3 2 YesDescription:
- ID: Identificator
- Cost: Annual cost of health insurance.
- Age: Age of a person.
- Gender (F, M).
- BMI: Body mass index.
- Kids: Number of children.
- Smoker: Does the person smoke? (No, Yes).
mydata$GenderF <- factor(mydata$Gender, levels = c("F", "M"))
mydata$SmokerF <- factor(mydata$Smoker, levels = c("No", "Yes"))library(dplyr)
mydata %>%
select(Cost, Age, BMI, Kids, GenderF, SmokerF) %>%
summary()## Cost Age BMI Kids
## Min. : 1141 Min. :18.00 Min. :17.80 Min. :0.00
## 1st Qu.: 5700 1st Qu.:31.00 1st Qu.:26.10 1st Qu.:0.00
## Median :10652 Median :43.00 Median :31.30 Median :1.00
## Mean :13789 Mean :42.34 Mean :31.15 Mean :1.25
## 3rd Qu.:14476 3rd Qu.:55.00 3rd Qu.:35.88 3rd Qu.:2.00
## Max. :60021 Max. :64.00 Max. :49.10 Max. :5.00
## GenderF SmokerF
## F:71 No :113
## M:69 Yes: 27
##
##
##
## library(GGally) # Load tools for scatterplot matrices
ggpairs(
mydata[c("Cost", "Age", "BMI", "Kids", "GenderF", "SmokerF")], # Variables to display
upper = list(
continuous = wrap("smooth", method = "lm", se = FALSE, alpha = 0.6) # Scatterplots with OLS line
),
lower = list(continuous = "cor"), # Pearson correlation coefficients
diag = list(continuous = "densityDiag") # Density plots on the diagonal
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 7),
axis.text.y = element_text(size = 7)
)
fit <- lm(Cost ~ Age + BMI + Kids + GenderF+ SmokerF,
data = mydata)
vif(fit) ## Age BMI Kids GenderF SmokerF
## 1.028518 1.024215 1.005644 1.004913 1.011513mean(vif(fit))## [1] 1.014961mydata$StdResid <- rstandard(fit)
hist(mydata$StdResid,
xlab = "Standardised residuals",
ylab = "Frequency",
main = "Histogram of standardised residuals")
shapiro.test(mydata$StdResid)##
## Shapiro-Wilk normality test
##
## data: mydata$StdResid
## W = 0.85838, p-value = 2.893e-10library(dplyr)
# Calculate Cook's distance
mydata$CooksD <- cooks.distance(fit)
# Define Cook's distance threshold
cd_threshold <- 4/nrow(mydata)
# Display apartments with Cook's distance above the threshold
influential_units <- mydata %>%
filter(CooksD > cd_threshold) %>%
arrange(desc(CooksD)) %>%
select(ID, Cost, Age, BMI, Kids, GenderF, SmokerF, StdResid, CooksD)
head(influential_units, 6) %>%
mutate(StdResid = round(StdResid, 3),
CooksD = round(CooksD, 3))## ID Cost Age BMI Kids GenderF SmokerF StdResid CooksD
## 1 1028 21595 23 18.7 0 M No 3.726 0.149
## 2 1231 60021 52 34.5 3 M Yes 3.624 0.146
## 3 1013 36580 61 33.3 4 F No 3.265 0.103
## 4 527 24060 19 30.6 2 F No 3.384 0.077
## 5 63 30167 64 24.7 1 M No 2.899 0.065
## 6 267 17180 40 19.8 1 M Yes -1.979 0.046mydata$StdFitted <- scale(fit$fitted.values)
library(car)
scatterplot(y = mydata$StdResid, x = mydata$StdFitted,
ylab = "Standardised residuals",
xlab = "Standardised fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riversols_test_breusch_pagan(fit)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## --------------------------------
## Response : Cost
## Variables: fitted values of Cost
##
## Test Summary
## -------------------------------
## DF = 1
## Chi2 = 31.01132
## Prob > Chi2 = 2.565275e-08library(estimatr)
fit2 <- lm_robust(Cost ~ Age + BMI + Kids + GenderF+ SmokerF,
se_type = "HC1",
data = mydata)
summary(fit2)##
## Call:
## lm_robust(formula = Cost ~ Age + BMI + Kids + GenderF + SmokerF,
## data = mydata, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) -12382.1 3763.39 -3.2901 1.280e-03 -19825.428 -4938.8
## Age 264.8 40.79 6.4916 1.521e-09 184.104 345.4
## BMI 304.1 98.36 3.0917 2.422e-03 109.557 498.6
## Kids 734.5 369.45 1.9882 4.883e-02 3.826 1465.2
## GenderFM -268.1 1053.46 -0.2545 7.995e-01 -2351.650 1815.5
## SmokerFYes 24380.5 1959.58 12.4417 4.233e-24 20504.825 28256.2
## DF
## (Intercept) 134
## Age 134
## BMI 134
## Kids 134
## GenderFM 134
## SmokerFYes 134
##
## Multiple R-squared: 0.7395 , Adjusted R-squared: 0.7298
## F-statistic: 37.91 on 5 and 134 DF, p-value: < 2.2e-16